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Investigation of cellular and network dynamics in the brain by means of modeling 
and simulation has evolved into a highly interdisciplinary field, that uses sophisticated 
modeling and simulation approaches to understand distinct areas of brain function. 
Depending on the underlying complexity, these models vary in their level of detail, in 
order to cope with the attached computational cost. Hence for large network simulations, 
single neurons are typically reduced to time-dependent signal processors, dismissing the 
spatial aspect of each cell. For single cell or networks with relatively small numbers of 
neurons, general purpose simulators allow for space and time-dependent simulations of 
electrical signal processing, based on the cable equation theory. An emerging field in 
Computational Neuroscience encompasses a new level of detail by incorporating the full 
three-dimensional morphology of cells and organelles into three-dimensional, space and 
time-dependent, simulations. While every approach has its advantages and limitations, 
such as computational cost, integrated and methods-spanning simulation approaches, 
depending on the network size could establish new ways to investigate the brain. In this 
paper we present a hybrid simulation approach, that makes use of reduced ID-models 
using e.g., the NEURON simulator — which couples to fully resolved models for simulating 
cellular and sub-cellular dynamics, including the detailed three-dimensional morphology 
of neurons and organelles. In order to couple ID- and 3D-simulations, we present 
a geometry-, membrane potential- and intracellular concentration mapping framework, 
with which graph- based morphologies, e.g., in the swc- or hoc-format, are mapped 
to full surface and volume representations of the neuron and computational data from 
ID-simulations can be used as boundary conditions for full 3D simulations and vice 
versa. Thus, established models and data, based on general purpose ID-simulators, can 
be directly coupled to the emerging field of fully resolved, highly detailed 3D-modeling 
approaches. We present the developed general framework for 1 D/3D hybrid modeling 
and apply it to investigate electhcally active neurons and their intracellular spatio-temporal 
calcium dynamics. 



Keywords: calcium dynamics, electrical stimulation, hybrid, parallel, detailed modeling, intracellular signaling, 
neuron, PDEs 



1. INTRODUCTION 

The level of detail, that can be achieved with experimental tech- 
niques in Neuroscience is ever increasing. Intracellular organelles, 
protein positioning, and messenger dynamics can be recorded 
in spatial and temporal sequences, (Spacek and Harris, 1997; 
Takahashi et al, 1999; Arellano et al, 2007; Chen et al, 2008; 
Popov et al., 2011; Burette et al., 2012). It is obvious that processes 
in neurons and the brain take place in three-dimensional space 
and that the three-dimensionality is employed by the biologi- 
cal system to regulate and differentiate highly complex signaling 
pathways in the brain. To address the role of time-dependent, 
three-dimensional signal processing and investigate the role of 
morphology on signaling from a computational standpoint, an 
inevitable step is to model critical processes in the brain in three 
space dimensions and in time. 



State of the art tools in Computational Neuroscience reduce 
three-dimensional problems to one-dimensional ones. For 
instance, simulators for the electrical signaling in neurons, e.g., 
NEURON (Hines and Carnevale, 2003) or Genesis (Bower and 
Beeman, 1997), solve a one-dimensional numerical problem in 
space. While this has great advantages in many applications, fore- 
most the computational speed of the methods that allows the 
simulation of large network activity, the drawback is the loss of 
modeling the intra- and extra- cellular space of neurons, and being 
able to include intracellular processes in a full three-dimensional 
resolution. The three-dimensional organization of neurons, e.g., 
the filigreed geometry of the endoplasmic reticulum, (Spacek 
and Harris, 1997) or the intra- and inter-ceUular organization of 
spines, (Murase and Schuman, 1999; Arellano et al., 2007; Chen 
et al, 2008; Tai et al, 2008; Popov et al, 2011), demonstrates the 



Frontiers in Neuroinformatics 



www.frontlersln.org 



July 2014 I Volume 8 | Article 68 | 1 



Grein et al. 



1D-3D hybrid modeling 



necessity for modeling intracellular processes in a highly detailed 
fashion in order to capture the underlying physical concepts of 
cellular signaling. 

The need for coupling different aspects of signal processing 
in neurons has been recognized in the past (Kerr et al, 2008; 
Wils and De Shutter, 2009; Andrews et al, 2010; Cannon et al, 
2010; Oliveira et al, 2010). In addition to projects focussing 
on the coupling of existing simulators for parallel computing 
architectures (Djurfeldt et al., 2010), the influence of spatial 
channel distribution on the electrical properties (Cannon et al., 
2010) or integrating reduced intra-cellular approximations of 
reaction-diffusion processes, (Resasco et al., 2012; Anwar et al, 
2013; McDougal et al., 2013a), we focus on the topic of how 
the three-dimensional intracellular architecture of neurons influ- 
ences intra-cellular signals and how the resulting models can be 
efficiently solved on different computing scales. Thus, the key 
factors are to incorporate an accurate morphology of the neu- 
ron including the three-dimensional intra-cellular architecture 
(which can include active and passive organelles), model a multi- 
ion system described by systems of partial and ordinary differen- 
tial equations (these can include the exchange mechanisms across 
plasma- and organelle membranes) and allow bi-directional cou- 
pling between one-dimensional membrane potential and three- 
dimensional intra-cellular signaling computations. The complex- 
ity of three-dimensional, detailed simulations of multi-ion sys- 
tems is efficiently handled by the simulation framework uG, 
where Finite Element or Finite Volume discretizations and multi- 
grid methods are used to solve the underlying large linear equa- 
tion systems on, potentially, massively parallel systems, (Heppner 
et al., 2013; Vogel et al., in press). Making use of uG, we propose 
a method to couple classical ID-simulators for the computation 
of membrane potentials, with detailed 3D-simulations within an 
on-line lD/3D-hybrid simulation framework. 

In this paper we demonstrate this newly developed approach 
using the example of intra-cellular calcium dynamics coupled 
to the membrane potential via voltage-gated calcium channels. 
While ID simulations can be carried out with classical simula- 
tors, such as NEURON or Genesis, the simulations in 3D are 
based on systems of partial and ordinary differential equations 
(PDEs and ODEs) describing distinct intracellular processes. In 
order to merge the two into a lD/3D-hybrid system, two things 
need to be done. In the first step one needs to generate from 
a ID morphology (a compartment model geometry) a surface 
and volume grid for numerical simulations in 3D and in the sec- 
ond step ensure a bi-directional coupling where the membrane 
potentials of the ID simulation are mapped onto the surface 
of the 3D problem as boundary conditions for the numerical 
problem and the computed intra-cellular ion concentrations are 
mapped back to compute updated membrane potentials. For this 
we developed automated tools to compute the necessary mor- 
phology representations of a neuron and to allow bi-directional 
coupling. 

For the neuron surface meshing in 3D, we employ a triangu- 
lation algorithm that makes use of the coordinates and radii of 
anatomical reconstructions, stored in e.g., hoc- or swc-format, to 
compute an equivalent surface mesh of the reconstructed neu- 
ron. We used TetGen (Si, 2009) to generate, from the surface grid. 



an intra- and extra-ceUular tetrahedral volume grid. In order to 
associate grid nodes from the triangular surface grid with nodes 
from the compartment model geometry we developed V,„2uG, a 
tool that uses a nearest neighbor-algorithm and kd-search algo- 
rithm to perform the mapping of membrane potential data onto 
the triangulated surface of the 3D-neuron. 

In this paper we chose a setup consisting of a multi- 
compartment model of a CAl stratum radiatum interneuron 
taken from Katona et al. (2011) and carried out simulations 
for electrical signal processing in NEURON. We then used the 
presented methods to couple ID and 3D models and simu- 
lated calcium dynamics in the full three-dimensional intracellular 
space under various parameter settings. Based on these examples, 
we introduce our methods and tools for lD/3D-hybrid modeling 
and show how existing models and data can be incorporated into 
highly detailed, three-dimensional simulations. 

2. RESULTS 

The methods described here couple one-dimensional electrical 
models and three-dimensional models for intra-cellular signaling. 
We chose the CAl interneuron from Katona et al. (2011) and its 
neuron morphology made available on NeuroMorpho.org (Ascoli, 
2006) and on ModelDB (MigHore et al, 2003). Simulations of 
the membrane potential dynamics in ID, i.e., on a compart- 
ment model level, were performed with NEURON (Hines and 
Carnevale, 2003; Carnevale and Hines, 2006), using a standard 
set-up defined in the Materials and Methods, section. Since 
intra-cellular processes are strongly regulated by calcium, e.g., 
(Milner et al., 1998; Bading, 1998; West et al, 2002; Clapham, 
2007; Tai et al, 2008), we chose calcium dynamics regulated by 
plasma membrane-located calcium channels with a given den- 
sity, thus modeling effectively a channel conductance density, and 
a diffusion-reaction process in the neuronal cytosol as a repre- 
sentative of three-dimensional, intracellular signaling in neurons. 
3D simulations were carried out in uG, Bastian et al. (1997); 
Vogel et al. (in press). Note, that this is a representative setup 
which is applicable to any other ID simulations and 3D intracel- 
lular processes. The coupling of both models here occurs on the 
level of calcium channels — these require the membrane poten- 
tial in space and time on the plasma membrane, the local intra- 
and extra-cellular calcium concentrations, as well as the geometry 
itself In this section we will introduce the models and the simula- 
tion set-up, methods for grid generation and membrane potential 
mapping, and wiU show simulation results for the described 
1D/3D hybrid simulation approach. 

2.1. THE 3D CALCIUM MODEL 

For this study we consider a calcium model on the continuum 
scale, including the following components: 

1. Morphology: The morphology and thus the computational 
domain is defined by a standard compartment model, e.g., in 
the hoc-format (see Supplemental Figure SI for an example). 
This morphology is then mapped to an equivalent three- 
dimensional computational domain. 

2. Membrane potential: The membrane potential is an input 
parameter for the calcium channel models and is computed 
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by the ID simulations and provided to the calcium channel 
models as input data. 

3. Calcium channels on the plasma membrane: Based on the 
models by Borg-Graham (Graham, 1999), we define N-/L- 
or T-type calcium channels (see Materials and Methods). 
Channel densities can be space-dependent, thus inhomoge- 
neous channel distribution is possible. 

4. Cytosolic calcium dynamics: In this study we consider diffu- 
sion of calcium and reaction of calcium with buffers in the 
cytosol. The computed calcium concentrations are mapped to 
the ID model to compute calcium-dependent currents. 

We can formulate the above points mathematically as an initial- 
value boundary problem for a diffusion-reaction model. To this 
end, let us denote the neuron geometry as f2, which is a compact 
subset of R^, such that 



• ^2 C with plasma membrane boundary F 

• f2 := f2 U af2 and 

• QndQ = 0 



9f2, 



which defines the space-time cylinder 

Z := Ux]0, r[c R^ X R+, V(x, t) e Z : T > t > 0. (1) 

Q is the problem-associated domain for the 3D model, recon- 
structed from the original compartmental model geometry 
defined, in our case, by a hoc file. Furthermore let c(x, t) be 
the calcium concentration function. Then the diffusion-reaction 
problem can be stated as 



dc(x, t) 
dt 



D- Ac{x, t) + R{c) inQ 



c(x, 0) = co(3c), in f2 
dc(x, t) 



dn 



0(V„,(5, f), X, f), on dU, 



(2) 
(3) 
(4) 



where D denotes the diffusion coefficient for cytosolic calcium, 
A := XiLi is the Laplace-Operator and R{c) a reaction 

term that depends on the calcium concentration c. Note that, 
the vector n denotes the direction perpendicular to the bound- 
ary surface. Equation (3) is the initial condition, i.e., a calcium 
distribution for t = 0 in the cytosol defined by the function cq (3c) . 
Function <I> defines a Neumann flux boundary condition, regu- 
lating the calcium flux through N-/L- or T-type calcium chan- 
nels, and thus also depends on the space- and time-dependent 
membrane potential V„, (x, f ). O in this study is specified by Borg- 
Graham model type calcium channels, which are listed in the 
Materials and Methods section. 

Cytosolic interaction of calcium with mobile and stationary 
buffers is governed by reaction equations of the type: 



B 



)n.B,„„l„-fe[Cfl^+] ■ [B^obile] — Kjf,B„Mic ' [CaB„ohile\ 

as well as the conservation law for the buffer concentration: 



[CaB„ 



ibilei 



B, 



mobile, total 



V^mohile\ 



(5) 



(6) 



The spatio-temporal dynamics of buffering molecules are 
described by a diffusion-process 



dB^otileix) 
dt 



Db,„ 



AB,nohileix) - B 



(7) 



For stationary buffers the diffusion coefficient DB^^f,^,^ can be 
set to zero. The parameter values of the model equations used 
throughout the paper are listed in Table 1. 

Note, that our framework is not limited to the above exam- 
ple of calcium/buffer dynamics. The employed multi-physics 
platform uG has been successfully used in a wide variety of appli- 
cations, ranging from, but not limited to, simulations of skin 
permeability in pharmaceutical applications (Hansen et al, 2008; 
Nagel et al, 2008, 2009; Muha et al, 2011) to groundwater flow 
studies (Grillo et al., 2010) and biogas reactor modeling (Muha 
et al, 2012). The latter demonstrates the use of the framework for 
large chemical reaction systems. The 3D intra-cellular model pre- 
sented here can thus be extended to include multiple ion species 
that are non-linearly coupled, resulting in a non-linear system of 
PDEs. 

2.2. COMPONENTS OF THE HYBRID FRAMEWORK 

As mentioned earlier, our hybrid 1D/3D framework consists of 
geometry matching and the coupling of ID membrane poten- 
tial simulations and 3D intra-cellular simulations. Computational 
grids for numerical simulations in 3D rely on a geometry-defining 
surface grid and a discrete representation of the computational 
domain, i.e., a volume grid. Approximation of unstructured 
domains is ideally done by a triangular surface and tetrahe- 
dral volume grid. These grids need to be generated using the 



Table 1 | Parameters used for the simulations using the 1D/3D hybrid 
frameworl<. 





NEURON 


uG 


dt [ms] 


0.1 


0.1 


stim.dur [msl 


10 




stim.amp [nA] 


0.1 




timesteps [#] 


10000 


10000 


Boundary Condition 




Vmfrom NEURON 


Calcium diffusion coefficient [\x.m^ ■ s"''] 


20-100 


20-100 


VGCC density Iixm"^] 


200-1000 


200-1000 


Buffer kon [iiM~^ ■ ms~^] 




0.09 


Buffer kaff [ms~^] 




0.24 


Buffer diffusion coefficient [(im^ ■ ms~^ 1 




0.043 


Total buffer concentration [|iM] 




20 


Extracellular calcium concentration [mM] 


1.6 


1.6 


Intracellular calcium concentration [[iM] 


0.1 


0.1 


Temperature [K] 


300 


300 



The calcium diffusion coefficient was varied across the ranges documented in 
Allbritton et ai. (1992) and the VGCC density was varied according to Pumplin 
etal. (1981); Eggermann et al. (2011). Values for the buffer kinetics and diffusivity 
were taken from Nagerl et al. (2000) as well as extracellular calcium concentra- 
tion from Egelman and Montague (1999) repsectively the intra-cellular calcium 
concentration from Helmchen et al. (1996); Maravall et al. (2000). 
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compartment geometry information, provided by the ID model 
as, e.g., a hoc- or swc-file. To provide an automated way for gener- 
ating surface and volume grids from anatomically reconstructed 
geometries, we developed novel tools and combined them with 
existing ones, which 

1. Integrate NEURON to perform ID simulations of membrane 
activity in neurons. 

2. Generate triangular surface grids from hoc-style neuron mor- 
phology files. 

3. Map ID simulation data of membrane potentials from 
compartment geometries onto the two-dimensional neuron 
membrane (VmluG) which are used as boundary condi- 
tions for membrane potential-dependent channels in the 
3D model. 

4. Integrate ion concentrations on the equivalent surface grid to 
feedback into the NEURON simulation. 

The general workflow of the 1D/3D coupling is depicted in 
Figure 1 . Using the graph-structure data contained in neuron 
compartment geometry files, we generate an equivalent 3D geom- 
etry with the compartment model graph as its "backbone." This 
procedure is fully automated and quality controlled and will be 
discussed in a forth-coming paper. Since our presented frame- 
work is not limited to the grid generator used here, we would 
like to acknowledge the efforts of other authors addressing sur- 
face grid generation from point-diameter information, (e.g., 
McDougaletal, 2013b). 



In Figure 2 we show the in hoc-style defined neuron over- 
layed by the corresponding three-dimensional neuron, which can 
be used for the fuU-spatial simulation of intra- cellular processes. 
The quality of the volume grid influences the numerical accuracy 
and stability. In particular the tetrahedral angles are a means of 
determining the grid quality (Deuflhard and Weiser, 2011). For 
the grid used in this paper we measured minimal and maximal 
aspect ratios of the the triangular surface grid ARt„ and the tetra- 
hedral volume grid ARtet respectively, yielding values between 
0.14 < ARtri < 0.86 and 0.21 < ARtet < 0.8. As a guideline one 
can state that aspect ratios above 0.1 result in grid elements that 
do not affect numerical stability (Shewchuk, 2002; Deuflhard and 
Weiser, 2011; Thompson et al, 2012). 

2.3. MEMBRANE POTENTIAL MAPPING— \//w2{/6 

In ID compartment models the membrane potential y,„ is com- 
puted in one node per cylindrical compartment. In the 3D setting 
each original cylinder is now represented by a segment of the 
triangular surface grid, thus the membrane potential from one 
cylinder node needs to be mapped onto a calculated cluster of 
nodes in the triangular surface grid. 

Associating each grid point y := iyi,y2,y3) of the 3D mor- 
phology with compartment model grid points x := {xi, X2, x^) 
(Figure 1) and the corresponding membrane potentials over 
the time-course of a simulation thus requires a mapping 

V„2uG: X R+ X R ^ X R+ X R (8) 
(x, t, VmJ (y, t, V„^) 



1 D Framework 

Set up 1 D multi-compartment model 
with NEURON: 

1 . Specify morphology with swc- or hoc-file 

2. Define electrical properties 

3. Specify intracellular quantities, e.g. calcium ions 
regulated by voltage-gated calcium channels 



3D Framework 

Setup 3D model with uG: 

1 . Create volumetric grid from underlying swc- or 
hoc-defined morphology 

2. Specify the biophysics of al! intra-cellular quantities and 
organelles (e.g. reaction-diffusion of calcium ions and mobile 
buffers with organelle interaction) 

3. Associate boundary nodes from volume grid with 
multi-compartment nodes 



feed W-setupinto 
hybrid model 



1 D/3D Hybrid Framework 



feed BD-setup into 
hybrid mode! 



ID Simulation - Compute new membrane potential with 3D-computed calcium concentrations 




- Compute voltage gated calcium fluxes with 1 D-computed membrane potential 
3D Simulation - Compute intracellular calcium dynamics (e.g. dependent on diffusion and 

reaction with mobile buffers, endoplasmic and mitochondrial calcium exchange) 



FIGURE 1 I Workflow of the 1D/3D hybrid framework. The workflow can 
be separated into 1 D and 3D modeling and simulation and into a set-up and 
simulation phase. In the set-up phase a general purpose simulator is used to 
define the one-dimensional problem (e.g., a hoc-geometry and set-up file). 
The 3D setup consists of generating a three-dimensional representation of 
the defined neuron and of specifying the cellular and intra-cellular 



components of the model problem. In the simulation phase 1 D membrane 
potential simulation data is computed and mapped to the 3D framework. 
There, the membrane potential data is included in simulating intracellular 
processes, such as calcium signaling. The computed calcium data is then 
mapped back to the 1 D problem, where it is used to compute 
calcium-dependent membrane fluxes. 
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B 



FIGURE 2 I Overlay of the multi-compartmental model defined by 
NEURON (red) and the equivalent tetrahedral volume mesh by uG 
(blue). (A) View of the whole neuron from Katona et al. (2011). Note, 
that close to the soma local optimization of the morphology leads to 
slight divergence between NEURON and uG morphology. Optimization is 
performed when the multi-compartment geometry model would cause 



unphysiological intersections of dendrites at the soma in 3D, due to the 
fact that all dendrites in the multi-compartment description share a 
common soma node. (B) Magnified views of ROIs color-coded by 
rectangular boxes within the overview for comparison, showing the 
compartment model geometry defining the backbone of the generated 
surface grid. 



that assigns a membrane potential V^y to every grid point y. V„,j, 
is set to the potential V^^ associated with the nearest neighbor 
3c := (xi, X2, X3) of grid point y with respect to all points con- 
tained in the hoc-file in every timestep t of the simulation. The 
distance measure for calculating the nearest neighbor can be any 
Minkowski metric and is interchangeable. In this study we used 
Algorithms 1 and 2 and chose the euclidean metric, or L2-norm, 
such that 

/d-l xl/2 

dist{p,q):=\^J2{p,-q,fj (9) 

Because the geometry information provided by the hoc-file may 
be very large, i.e., the number of provided points n is large, we 
implemented and tested two strategies for determining the near- 
est neighbor. The first is an exact solution of the problem by space 
partitioning with multi-dimensional binary search trees, k-d trees 
with a search complexity 0(n ■ login)) (see Bentley, 1975 as well 
as Figure 3 for an example), the second strategy is the exact solu- 
tion by pairwise comparison with a search complexity 0{n^). 
Note that we consider a three-dimensional example here, but the 



mapping process is dimension-independent and can be done for 
arbitrary dimensions by means of the underlying C-|~|- library for 
multi-dimensional (binary) search trees (Mount and Arya, 2012), 
a library that is used in Vm2uG. 

2.4. ALGORITHMS 

The pairwise comparison is an exhaustive search that requires no 
additional data structure but 0(m • m) iterations which is imprac- 
tical if the number n of provided points on the grid is large, and 
for the number of corresponding points in the hoc file holds n ~ 
m, thus arriving at quadratic runtime complexity O(m^). For large 
M, we therefore use the method of space partitioning (k-d trees. 
Figures 3B,C) which leads to an improved average runtime com- 
plexity for a query of log (^), with d being the space dimension. 
With this approach we reach super-linear runtimes after the ini- 
tial tree has been built. In contrast to the 0{n^) method, we need 
an additional structure, which is justified because of the overall 
speedup compared to the naive approach. The overall runtime is 
dominated by the utilized search algorithm used during the ini- 
tial build of the tree and the balancing of the tree itself, yielding a 
complexity Oin ■ log(«)) (Wald and Hvran, 2006; Gormen et al, 
2011). The worst case arises iin<^2'^ is satisfied and the runtime 
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FIGURE 3 I Basic operations. (A) General description of the nearest 
neighbor search task in 2D/3D: Let S be a set of #n points with dimension 
d := 2, 3. Let q be a query point. The task is to find a point r (green) in the 
set S of all points (black), which is closest to the query point q [center of 
sphere in (A) on the right] within a given bounding box, i.e., within a circle 
(black) or sphere. (A) Pairwise comparison of all points in a radius search. B: 
Space-partitioning approach using k-d trees (multi-dimensional binary 
search trees). (C) An exemplary multi-dimensional binary search tree of 
dimension d = 2 (enhanced and modified according to Wikipedia, 2012), 
representing the space partition in (B). 



Algorithm 1 | Linear interpolation 

q ^ query point 
2; NN +- getNearestNeighbor(q) 

iDist +- NN.getDistanceO 
4: if dist < cutoff then 

NNs *- getNearestNeighbors(k) 
6: for i = 1 ^ k. j = 1 ^ k A i j!= j do 

g *- constructLine(NNs[il, NNsljl) > create straight line 

8: q' -s- projectOnLine(q, g) > orthogonally project onto line 

if getDistance(q', q) < iDist then 
10: iDist ^ getDistance(q', q) 

Vm ^ g.interpLinear(q') > linear interpolation 

12: end if 

end for 
14: end If 



complexity increases per query to O (^d ■ n^^ s'^ (Lee, 1997). 

Yet, the curse of dimensionality is not an issue in the presented 
application, since our geometry is defined in three-dimensional 
space. 

V„2uG offers linear or bilinear membrane potential interpo- 
lation (Algorithms 1 and 2). This can be useful if the compart- 
mental representation of the neuron is very coarse and map- 
ping distances become large, i.e., 3D surface grid points lie far 
away from the computed nearest neighbor on the compartment 



Algorithm 2 | Bilinear interpolation 

q ^ query point 
2: NN ^ getNearestNeighbor(q) 

iDist ^ NN.getDistanceO 
4: if dist < cutoff then 

NNs ^ getNearestNeighbors(k) 
6: for / = 1 ^ *:,/= 1 k. 
l=^^kAiJ^jJ^k^l do 

det *- calculateDeterminant > calculate determinant 

8: p constructPlane(NN[il, NN[jl, 

NN[kl, NN[I1) > create plane 

q' ^ projectOnPlane(q, p) 
10: if getDistanoe(q', q) < iDist then > orthogonal project onto plane 

iDist ^ getDistance(q', q) 
12: Vm ^ p.interpBilinear(q') > bilinear interpolation 

end if 
14: end for 
end if 



geometry. Interpolation is then carried out over n e pseudo 
nearest neighbors to compute an optimal approximated value for 
the membrane potential which is then assigned to the given grid 
point. 

2.5. WORKFLOW AND BIDIRECTIONAL COUPLING 

We use a direct data coupling mechanism termed a type A problem 
in Heterogenous Multiscale Modeling, see (Weinan et al., 2003), 
utilizing an online algorithm within uG. The NEURON library is 
used as a plugin/shared library within the uG project. The sim- 
ulation workflow is defined in uG by a lua-script (lerusalimschy 
et al., 2013) (see Figure 4 for the script used in this paper). We 
provide fully flexible control over NEURON within uG by means 
of a custom API or by directly including the hoc-script (which is 
then executed by the HOC language interpreter). For the latter we 
developed an additional C-|— I- wrapper. 

The workflow defined in the lua-script performs the steps 
illustrated in Figures 1, 4. The numerical procedures are defined 
individually for NEURON and uG (see Materials and Methods), 
which is necessary since the model equations for the ID and 3D 
model are of different types. The only coupling requirement is 
that both simulations can synchronize their data at specific time 
points, which is guaranteed by our framework. 

Data exchange between ID and 3D model is bidirectional, 
where membrane potential data is passed from ID to 3D via the 
mapping algorithm presented in the previous section. The com- 
puted ion concentrations in the 3D model are then passed from 
3D to ID in the following way: 

Consider a cylindrical compartment C, of the ID model. By 
9(Q n 9^2) we denote the portion of the cell surface d^l that is 
associated with compartment C,. We can then compute the intra- 
cellular calcium concentration for C, as 

2+ ll9C,il f 

\Ca Ir = / Urn2+ix) 

'''' ||3(Q 0 3^)11 Ja(c,n3a) 

forie 9(C, n 9^2) and Vf = 0,...,N (10) 



Frontiers in Neuroinformatics 



www.frontiersin.org 



July 2014 I Volume 8 | Article 68 | 6 



Grein et al. 



1 D-3D hybrid modeling 



1 BliJ 

2 -- NEURON setup 

3 --]] 

4 -- Command line arguments 

5 hoc_stim = util.GetParam["-hoc_stiin", "~/hoc/stim.hoc") 

6 hoc_geom = util.GetParam("-hoc_geoin", "~/hoc/geom.hoc") 

7 hoc_dt = util.GetParamNumber("-hoc_dt", 0.1) 

8 hoc_tstop = util . GetParamNumberC" -hoc_tstop" , 1.0) 

9 hoc_tstart = util . GetParamNumberC -hoc_tstart" , 0.0) 

10 hoc_finit = util . GetParamNumberC -hoc_f initialize" , -75.0) 

11 -- Initialize 

12 HyHocSetup = HocSetupO 

1 3 HyHocSetup : load_geom (hoc_geom) 

14 HyHocSetup : load_stim (hoc_stim) 

15 HyHocSetup ; setup_hoc (0, hoc_tstop, hoc_dt, hoc_finit) 

16 -- Summary 

17 HyHocSetup;print_setup(true) 

18 --[[ 

19 -- UG setup 

20 --]] 

21 -- Variables 

22 bg_conductivity = IQOO.O * 0.06 

23 bg_density = lOQO.O 

24 bg_type = "N" 

25 delta_t = le-4 

26 numTimeSteps = 10000 

27 -- Channel setup 

28 function VGCC_Channels (Ca_i, Ca_o, x, y, z, t) 

29 local b = BorgGrahamChannels () 

30 b ; install_channel_gates ( bg_conductivity, bg_density, bg_type) 

31 if (t == 0) then 

32 b;calc_current_at_start(-75) 

33 else 

34 b;get_current(Ca_i, Ca_o, y, z, t) 

35 end 

36 return b;get_flux() 

37 end 

38 -- Chemical model specification 

39 ElementDiscretization = ConvectionDiffusion ("ca" , "Inner", "fvl") 

40 ElementDiscretization ; set_dif fusion (10. 0) 

41 ElementDiscretizationObstacle = ConvectionDiffusion ("ca" , "Obstacle", "fv") 

42 ElementDiscretizationObstacle ; set_dif fusion [1.0) 

43 neumannDisc = NeumannBoundary ("ca" ) 

44 neumannDisc ; add (VGCC_Channels, "Input", "Inner") 

45 u = GridFunction (ElementDiscretization, ElementDiscretizationObstacle) 

46 HyHocSetup : adiust_resistivities_obstacle (u, "ca", "Input") 

47 HyHocSetup ; calculate_and_set_initial_f eedback (u, "ca", "Inner") 

48 --[[ 

49 -- Run simulation 

50 --]] 

51 for step = 1, NumTimeSteps do 

52 PDESolver:apply(u) 

53 HyHocSetup ;calculate_and_set_feedbacks(u, "ca", "Inner") 

54 PDESolver; save_solution_vtk (u, file) 

55 end 



FIGURE 4 I Script file illustrating the setup of the NEURON and uG 
model, here divided into a NEURON setup phase, a uG setup phase, 
chemical model specifications and the simulation run control. One can 

make use of the hoc interpreter features within the lua script files that define 



the uG-setup. Note that one can specify the geometry and stimulation 
protocol in distinct files and combine this with using the hoc interpreter from 
uG to refine the setup by statements of the form: 
HocSetup:execute("command string "). 



where U(j^2+{x) is the calcium concentration in node x of the 
3D grid and N is the number of compartments representing 
the ID geometry. ||.|| symbohzes the size of 9C, and 9(Q fl 9^^) 
respectively. 

The factor accounts for the fact, that the surface size 

||()(Cin3S2)|| 

of 9(Cj n 9f2) and the surface of Q are not necessarily identical. 
Computing the calcium fluxes in uG and in NEURON using 



the values computed by Equation (10) shows good agreement 
between the two models, see Supplemental Figure S2 as well as 
Supplemental Figure S3. 

2.6. SIMULATIONS USING THE 1D/3D HYBRID METHOD 

We applied the described method for coupling ID and 3D sim- 
ulations to a CAl interneuron from Katona et al. (2011), taken 
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FIGURE 5 I Geometry specifications. (A) Three-dimensional neuron from 
Katona et al. (2011) used in the simulations. The box highlights the region of 
interest used in (B,C). (B) Triangular surface grid of a dendritic segment 
including an intracellular obstacle (yellow). The arrow in (A) indicates the 
position of the obstacle. (C) Triangular surface grid of a dendritic segment 
without the intracellular obstacle. The surface grid can remain unaffected 
by the insertion of an intra-cellular object. 



Table 2 | Scaling. 



Problem size [# nodes] 


Procs [#] 


Runtime [s] 


4000 


1 


59 


16,000 


4 


58 


64,000 


16 


47 


128,000 


32 


49 


256,000 


64 


49 



The W/3D hybrid framework exhibits weal< scaling on the in-house cluster {32 
nodes with each providing an Intel Core 17 CPU @ 2.67 GHz running Linux). 
We solved the presented model equations using a Newton solver for the non- 
linear part and used a BiCGStab method as iterative solver combined with 
ILU-smoothened geometric multigrid preconditionen Each runtime is the aver- 
age of n= 10 runs each, where each run was a simulation of 1 /10 s real time and 
a numerical time step width of 7 ms. 

from the database ModelDB (Migliore et al, 2003). We set up 
a NEURON simulation with stimulation protocol (see Materials 
and Methods) and generated a three-dimensional geometry for 
calcium simulations that use the membrane potential data from 
the NEURON simulation (Figure 5) and feeds back the com- 
puted concentrations (see previous section). An overlay of hoc- 
and ugx- (the uG grid format) geometry shows good agreement 
between the two. 

The scientific gain from 3D simulations lies in a physically 
detailed representation of morphology and the underlying bio- 
physical processes in neurons. These processes are regulated by 
space and time-dependent variations of parameters like chan- 
nel densities, cytosolic and organelle architecture, diffusivity and 
of course cellular and organelle morphology. Previous work 
has demonstrated the importance of investigating the spatial 



organization of ion channels and their stochastic behavior, 
(Oliveira et al., 2010; Brandi et al, 2011), in part using a multi- 
scale modeling approach. Brandi et al. (2011) for instance present 
in their abstract a platform for combining NeuroRD (Blackwell 
et al, 2013) and MOOSE (Ray and Bhalla, 2008). The frame- 
work employed here can be viewed as a natural extension of 
previously published multi-scale approaches, with a strong focus 
on resolving the intra-cellular space with additional accuracy, 
rather than interpreting the intra-cellular space as a homoge- 
nous space for reaction-diffusion processes. Obstacles or active 
organelles, such as calcium-relevant stores (e.g., the endoplas- 
mic reticulum or mitochondria), can be easily integrated into 
our hybrid framework. In addition to this, our framework aims 
at simulating entire neurons up to small networks with a high 
level of cellular and intra-ceUular detail, rather than investigat- 
ing e.g., the single molecule level by Smoldyn Andrews et al. 
(2010). To accomplish this, we model at the continuum scale 
and make use of the massively scalable code uG, (Heppner et al., 
2013), demonstrated to perform ideal weak scaling up to 64k 
processes on the Hermit Supercomputer. In Table 2 we show 
ideal weak scaling for the presented simulations on our in-house 
cluster (see Materials and Methods), solving a 3D problem with 
256,000 grid points in under 1 min on 64 processors. In com- 
parison to other simulators, such as VCell or MCeU, the VCell 
website lists a 2D benchmark with a problem size of 5000-10,000 
grid points and a PDE-system with 4—5 PDEs being solved in 
roughly 20 min. In Balls and Baden (2004), MCell scaling stud- 
ies up to 64 cores show a scaling efficiency of 85-92%, solving a 
benchmark problem also in roughly 20 min. Since the underlying 
models of the compared simulators differ, a direct comparison 
is not necessarily valid. The presented hybrid method, based on 
a continuum model is applicable ideally for whole- and multi- 
cell simulations, while e.g., MCeU uses particle-based methods 
that so far do not go beyond spine simulations, i.e., microdomain 
simulations. 

In order to demonstrate some advantages of a 1D/3D hybrid 
approach, we defined a standard 1D/3D simulation protocol (see 
Materials and Methods) and then varied parameters, such as cal- 
cium diffusivity (Allbritton et al., 1992), voltage-gated calcium 
channel (VGCC) densities (Pumplin et al., 1981; Eggermann 
et al., 2011), and added intra-cellular obstacles to the cytosol. 

Figure 6 and the supplemental movies show calcium simu- 
lations in a region of interest using the 1D/3D hybrid method, 
with the calcium diffusion set to 100 jim^/s, a VGCC-density of 
1000 |im^^ and no obstacle present in the cytosol (see Figure 5 
for the neuron geometry used in the simulations). This exam- 
ple shows how the NEURON-calculated membrane potentials 
are mapped to the 3D model, where they regulate VGCCs 
and detailed calcium dynamics can be simulated in fuU three- 
dimensional space. The calcium concentrations are then mapped 
back to NEURON according to the previous section. 

2.6.1. Varying the diffusion coefficient 

We then varied the diffusion coefficient for intracellular calcium 
between 10 and 100|i,m^/s, which is within the experimen- 
tally observed ranges, (Allbritton et al, 1992), while the VGCC 
density was set to 1000 |j.m^^ (Figure 7A). Figure 7B shows a 
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H ^ 





FIGURE 6 I Time series for a simulation of intra-cellular calcium signaling 
with a calcium diffusion coefficient set to 100 M-m^/s, VGCC density of 
1000 (i,m~^ and no obstacle present in the dendrite. The basal calcium 
concentration is set to 100 nM. (A-F) show the propagation of calcium in a 
dendrite, triggered by the the stimulation protocol listed in Table 1. 
Membrane potentials were computed in NEURON and mapped to the plasma 
membrane in order to compute the VGCC-dynamics regulating the calcium 
exchange. Calcium concentrations were then mapped back to compute the 



Time (t) [s] 
Location -pm between -center 

calcium-dependent membrane potential in NEURON. (A-F) covers Is of rea 
time. Time points shown are at 0, 0.01, 0.02, 0.03, 0.1, and 1 s. Note, the 
maximum single-channel conductance was set to 60 pS according to Graham 
(1999). (G) Highlighted measuring points close to the plasma membrane (pm) 
in between and at the center used in (H). (H) Evolvement of the calcium 
concentrations at three different points over a period of 20 ms. Note that the 
VGCCs close at the peak amplitude and the calcium profile quickly adjust to a 
uniform value based on the intra-cellular diffusion. 



linear dependence to a global variation of the cytosolic calcium 
diffusion coefficient. Note, that the diffusion of calcium does not 
have to be isotropic, but can have different diffusion properties in 
the three space-dimensions. The diffusion coefficient then needs 
to be formulated as a diffusion tensor 



D := 



'40 0 
0 dy 0 

,004 



which can be defined in the uG- workflow (Vogel et al., in press) 
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Time (t) [s] 
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VGCC density [urn"'] 
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Time (t) [s] 
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Time (t) [s] 
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FIGURE 7 I Parameter variation. (A) Variation of tine diffusion 
coefficient regulates the intracellular calcium profile (using a constant 
VGCC density). The concentration profiles computed by the hybrid 
framework are validated in Figure S3 using an alternative, reductionist 
approach. (B) Quantification of the peak amplitudes in (A), = 0.9461 
and p = 0.003, showing a linear dependency between the diffusion 
coefficient and the peak amplitude. (C) Density variation of VGCCs 
sharpens the concentration profile in the considered dendrite (using a 
constant diffusion coefficient). (D) Quantification of the peak amplitudes 
in (C), = 0.99 and p= 0.018 showing a linear dependency between 
VGCC density and peak amplitudes. Calcium profiles are validated by 
calcium imaging data, Helmchen et al. (1996); Kits et al. (1997); Maravall 
et al. (2000). VGCC densities are varied according to documented 
physiological intervals up to 10000/(im^, as reported in Eggermann 
et al. (2011). (E,F) Placing obstacles of various sizes in the dendrite, 
from 0 to 100 % (see Table 3) affects the locally available calcium 
concentration, (E) before the obstacle and (F) behind the obstacle. 



2.6.2. Varying channel densities 

Next, we decided to vary the VGCC densities on the plasma mem- 
brane. We placed N-type VGCCs (see Materials and Methods) 
homogeneously along the dendrites at densities between 200 
and 1000 |xni^^ (Figure 7C). Due to an increased channel 
density we observe changes in the rise and decay times, as 
well as different peak amplitudes in the calcium profile. In 
Figure 7D we see that the peak amplitude of the intra-cellular 
calcium signal depends linearly on the VGCC density. Global 
and time-dependent changes in the VGCC distribution of 



the neuron can thus finely regulate the intra-cellular calcium 
code. 

2.6.3. Including intra-cellular obstacles 

Where in a ID model one needs to average the intra-cellular 
space and approximate morphologies, the 3D approach allows a 
full and detailed investigation. To demonstrate how intra-ceUular 
objects can be included into a simulation, we investigated calcium 
dynamics in the observed dendrite with and without an intra- 
cellular obstacle. For this we placed an obstacle (this could be 
e.g., a mitochondrion or endoplasmic reticulum) in the center of 
a dendrite (see Figure 5), defined as a cylinder that is 0.14 |xm 
long and has a variable circumference Cobstade-> embedded in a 
dendrite of average local circumference Cnemite- We then define 
the hindrance factor H of the intracellular obstacle with respect to 
the surrounding cylinder for that particular section as: 

H:= e [0,1] (11) 

^cylinder 

The obstacle sizes used in our simulations are listed in Table 3, 
which correspond to a dendritic occupancy between 0 and 100% 
(100% corresponds to a full blockage in the the dendrite). We 
chose the obstacle location to be right before the most distal 
bifurcation point of a certain neurite which emerges — and can 
be visually traced in Figure 5A — from the left most bottom of 
the soma itself Note, that the obstacle is placed in the middle 
of the considered neurite in the present use case. Intra-cellular 
organelles, occupying parts of the cytosol, affect the intra-cellular 
resistivity. To account for this in our ID model, one can either 
directly modify the axial resistance by changing the intra-cellular 
resistance or by changing the diameter of the cylinders that har- 
bor an organelle. The axial resistance Ra in a cylinder with length 
/ and diameter d is computed by: 

4' , . 

Ra = Rr--pr- (12) 

where R, is the intra-cellular resistance. The effective diameter can 
be calculated for all affected compartments by: 

^new — ^old ^obstacle (13) 

Our framework automatically checks which compartments are 
occupied by organelles and adjusts the diameters accordingly. We 
observed that local intracellular signaling is affected by obsta- 
cle size and position (Figures 7E,F). In more complex situations, 
e.g., space occupancy by mitochondria, endoplasmic reticulum or 
spine apparatuses in dendrites or spine necks, the inclusion of 
these functional organelles in a three-dimensional model might 
be necessary for explaining the intracellular dynamics of the neu- 
ron. Active intra-cellular organelles can be easily added within the 
presented framework. The user needs to define the geometry of 
the organelle (which can be a detailed three-dimensional recon- 
struction) and the biological exchange mechanisms across the 
organelle membrane (e.g., IP3-receptors, ryanodine receptors, 
SERCA pumps, sodium-calcium exchangers etc.). 
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Table 3 | Hindrance factors. 



Hindrance [%] 


Obstacle volume [|<.m'^] 


Obstacle surface area [).im^] 


0 


0.0000 




0.0000 


25 


0.0022 




0.3040 


50 


0.0344 




1.9000 


75 


0.0945 




3.7240 


100 


0.4076 




8.6401 



Listed above are the surface and volume sizes of tlie obstacles that were placed 
In the dendrite for the simulations shown In Figure 7 with corresponding hin- 
drance factors. A hindrance factor of 100 % represents a full blockage of the 
dendritic cross-section. 



3. DISCUSSION 

Modeling and simulating cellular and network processes are 
constrained by the complexity of the computational problem 
and the availability of experimental data to validate the mod- 
els. Experimental techniques are evolving, shedding new light on 
the filigreed functioning of neurons and networks. Spatial and 
time resolutions are becoming fine enough to resolve positions 
of single proteins and receptors at synapses and the intra-cellular 
domain can be investigated at increasing levels of detail. On the 
other hand computational complexity increases when a model 
includes biological and spatial detail, resulting in a highly detailed 
three-dimensional model of neurons. Computational complexity 
further increases when not only one neuron, but entire networks 
need to be modeled. 

To cope with network complexity, several approximation 
methods for cellular and network function have been devel- 
oped over the last decades. These approximations yield zero- 
or one-dimensional models in space, describing neurons with 
point- or compartment neuron models. Numerous specialized 
and general purpose neuron simulators incorporate these basic 
abstraction techniques (e.g.. Bower and Beeman, 1997; Hines 
and Carnevale, 2003; Gewaltig and Diesmann, 2007) and allow 
the user to simulate large and complex neural network struc- 
tures. These models and simulators were developed and evolved 
not only around the limitations of computational power, but 
also around the availability of experimental data to validate the 
models. The evolution of experimental methods, computational 
resources and numerical mathematics motivates the develop- 
ment of neuron models that include greater physical detail. 
This brings us to three-dimensional models that include the 
detailed morphology of neurons, either by modeling physical 
processes using stochastic (e.g., Andrews et al., 2010; Oliveira 
et al., 2010) or deterministic models. We make use of the lat- 
ter approach with the additional feature of being able to inte- 
grate obstacles and intracellular organelles, which yields models 
formulated on the continuum scale by means of systems of par- 
tial differential equations. Solving these equations numerically 
requires advanced mathematical tools. One general purpose plat- 
form is uG (Bastian et al. 1997), which we used in this paper. 
uG offers numerical discretization methods and efficient solvers 
for systems of partial differential equations on highly unstruc- 
tured grids and runs on massively parallel systems (Heppner 
et al., 2013; Vogel et al, in press). As shown in the past, 
these numerical tools are applicable and efficient on a cellular 



level or in a networks with relatively small numbers of neu- 
rons (Xylouris et al., 2007; Wittmann et al., 2009; Xylouris, 
2013). 

Given the limitations in computing resources, it is currently 
not feasible to model and simulate large networks of neurons 
with full single cell, three-dimensional detail. To make use of 
the advantages of highly detailed three-dimensional models and 
to cope with the complexity that comes with modeling large 
networks of neurons, we developed a method to couple state 
of the art general purpose neuron simulators, e.g., NEURON, 
with three-dimensional models of single neurons. This 1D/3D 
hybrid method includes automated tools to either reconstruct 
three-dimensional neuron morphologies from raw microscopy 
data (Broser et al, 2004; Queisser et al, 2008; Jungblut et al, 
2011), from anatomically recorded data (Wolf et al., 2013) or 
from graph-structure morphologies as used, e.g., in the NEURON 
Simulator in the form of hoc-files or swc-files. 

In this paper we introduced a framework for geometry and 
membrane potential and intra-cellular ion concentration map- 
ping between ID simulations and the equivalent 3D model. For 
this we used the NEURON simulator to compute electrical sig- 
nals in a compartment neuron and uG to simulate intracellular 
calcium dynamics in 3D. With this approach we were able to 
exploit the modeling and computational advantages that a gen- 
eral purpose simulator for large networks brings, as well as the 
necessary tools to investigate intra-cellular (or even extra-cellular) 
processes on a very fine scale. By including the detailed mor- 
phology of neurons — which can be subject to temporal adaption 
due to neuronal activity (Silver et al., 1990; Korkotian and Segal, 
1999; Muller et al, 2002; Van Aelst and Cline, 2004; Hayashi and 
Majewska, 2005; Tada and Sheng, 2006; Wittmann et al, 2009; 
Kanamori, 2013) — the interplay between cellular/organelle mor- 
phology and cellular function can be systematically investigated. 
Using defined stimulation protocols in NEURON, we ran ID- 
simulations mapping the results onto the three-dimensional mor- 
phology as boundary conditions for a cellular calcium model 
and vice versa. The calcium model consisted of a local distribu- 
tion of N-type voltage-gated calcium channels, regulated by the 
membrane potential and intracellular diffusion and reaction of 
calcium ions with a mobile buffer. After presenting the function- 
ality of the 1D/3D hybrid method on a reference parameter set, 
we varied the density of voltage-gated calcium channels, the cal- 
cium diffusion coefficient and introduced intracellular obstacles. 
The results show that for one, the presented method is designed in 
a general fashion and is thus applicable to a broad range of neuro- 
biological questions and that the effect of intra-cellular obstacles, 
locally changing channel densities and cytosolic diffusivity have a 
substantial effect on intracellular signals. 

For future research, problem-specific models can be used 
within the presented framework. For instance, not only single 
neurons, but entire networks could be simulated on the ID level, 
where a small subset of "key neurons" in the network can be 
resolved in fuU three-dimensional detail (Xylouris et al., 2007). 
Furthermore intra-cellular models for different ionic species, 
detailed models of intra-ceUular organelles, such as mitochondria 
and endoplasmic reticulum calcium exchangers, protein synthesis 
and synapses could be included on the 3D scale. 
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4. MATERIALS AND METHODS 
4.1. THE 1D ELECTRICAL MODEL 

To simulate the electrical activity of the model neuron, we used 
the classical Hodgkin-Huxley equations (Hodgkin and Huxley, 
1952) including sodium, potassium, calcium and a leakage cur- 
rent: 



dV 



of 

+gL(V - El) + Ica, 



(14) 



where Vm is the membrane potential, Cm the membrane capac- 
ity, gMa> gK and gL are the maximum single channel conductances 
for sodium and potassium channels, as well as a leakage current. 
Ek, Enb and El are the reversal potentials. Ica is the calcium cur- 
rent specified by Equation (26). The time and voltage-dependent 
gating variables n,m,h determine the gating behavior of the 
channels and are computed by the ordinary differential equation 



dx Xo 
9t ~ ~ 



(15) 



where x = n, m, h. The parameters x^c and in the equations 
above are computed by 



and 



1 



Olx + fix 



(16) 



(17) 



where Tx,o = 0 for x = n, m, h. The values for Ux and fix are 
derived from (Hodgkin and Huxley, 1952): 



0(n 


omiv+w) 

~ V-IO 

£10 — 1 


(18) 




V 

£80 

~ 8 




Pn 


(19) 


dm 


0.1(y-F25) 

~ V + 25 

£10 — 1 


(20) 


Pm 


V 

= 4ei8 


(21) 




= 0.07e2B 


(22) 


Ph 


1 


(23) 
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Solving a multi-compartment model where for each cylindrical 
compartment an equation of the type 



dV„ 
dt 



~f~ ^electrode 



(24) 



is computed using the NEURON simulator. Numerically, the aris- 
ing sets of ordinary differential equations are solved with the 



Crank-Nicolson method implemented in NEURON. We chose 
the CAl stratum radiatum hippocampal interneuron published 
in Katona et al. (201 1 ) and ran the following stimulation protocol: 
We inserted a current clamp at the mid of the soma and stim- 
ulated with amplitude of 0.1 nA with timestep dt = OA ms for a 
total time of 10 ms. We then simulated 1 s of membrane potential 
propagation in NEURON evoked due to that stimulation. 

4.2. VOLTAGE-GATED CALCIUM CHANNEL MODELS 

In order to include voltage-gated calcium channels, we used the 
Borg-Graham model described for different ion channel types in 
detail in Graham (1999). In our setup we included N-type cal- 
cium channels gates, though the implementation allows us to also 
include other types (e.g., L-/T-type gates). We define a mapping 
of the calcium fluxes calculated by the VGCC-model onto the 3D 
morphology: 



Borg Graham : 



. : {xb, t, v„.) _F (25) 



The dynamics of calcium ionic fluxes are described in the fol- 
lowing way (Graham, 1999) and are computed inside the ID- 
simulation loop: 

IcaiV, t, [Ca]„ [CaW = G(V, t)FiV, [Ca]„ [Cfl]„) (26) 

G describes the properties of the gating functions, in particular 
the open and close state-probabilities of the channels and is spec- 
ified in Equation (28). F can be computed using the GHK model, 
yielding 



F{V, [Ca]„ [CaW 



PCa- 



RT 



[Ca- 



2+n 



[C. 



,2+1 



exp( 



1 — exp( 



(27) 



where pca is the permeability and for N-type calcium channels is 
set to the value 10^^ err? js. F denotes the Faraday constant, R the 
gas constant and T the temperature. For the VGCCs used in the 
simulations, the gating function G is defined as 



G(y, t) = fc(y, t)i^{v, t) 



(28) 



Here, k and / fulfill the ordinary differential equation (15) and eqs. 
(16, 17), using the following values: 



UxiV) = Kxsxp 



ZxVxiV- Vl/2.x)F 

RT 



Px(V) = Kx exp ( — I (29) 



The parameters V"i/2. z, K, K, tq also depend on the particular 
channel and are documented in Graham (1999). We set the 
values accordingly: 
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Vi/2,i := - 21 mV 
Zk :=2 
Yk :=0 
Kj^ :=1.7 ms 
Xk,Qi :=1.7 ms 
Vxi2,i ■■= - 40 mV 
zi :=1 

Ki :=70 ms 
T; 0 :=70 ms 

It is sufficient to solve Equation (15) using the first order explicit 
Euler method, in which the required membrane potential is taken 
from the membrane potential mapping (V„,2uG) in each time 
step. 

4.3. NUMERICAL APPROACH 

The calcium model used in the presented simulations is based 
on a cytosolic diffusion equation with boundary conditions that 
include the plasma membrane calcium exchange mechanisms, 
i.e., Neumann boundary conditions representing an ionic flux 
over a given surface area of the plasma membrane. This can be 
stated as: 

9 - , 
V- (DVm(x, f)) = — m(x, f)inf2 cM (30) 
dt 

with boundary conditions 
du 

— = g(Channels) (31) 
dn 

where u denotes the intracellular calcium concentration, D the 
diffusion tensor (a 3 x 3-Matrrx) and g{Channels) the calcium 
flux depending on the channels defined in the simulation (see 
previous section). To solve Equation 30, we discretized in space 
by means of first order finite volumes, treating the time derivative 
by an implicit Euler method. 

For this we discretized the plasma membrane geometry with a 
triangular surface grid and the volume of the cell with a tetrahe- 
dral volume grid. In the following we denote the discrete volume 
of the neuron as Q. and each arising control volume box as B, , i = 
!,...,«. Therefore Q = [J" Bi. The Finite Volume method (cf 
Eymard et al., 2000) defines an ansatz space V/, C H^{^) which 
is spanned by piecewise linear ansatz functions A, (3c,) = Sy and 
is used as the approximation for the solution u. The approxima- 
tion of u is then defined as the linear combination of these ansatz 
functions, resulting in a system of linear equations 

n 

J^AgUj = bi (32) 
j 

where A is the stiffness matrbc. To solve the linear system of equa- 
tions, we applied a geometric multigrid solver (smoother ILU) as 



a preconditioner along with BiCGstab as the base solver for the 
linear part (inner loop) and a Newton solver for the non-linear 
part (outer loop), see e.g., (Hackbusch, 1985, 2003). The numeri- 
cal setup and simulations were carried out in uG [see Bastian et al. 
(1997); Vogel et al. (in press)]. 

ACKNOWLEDGMENTS 

Numerous gratitude is given to Andreas Vogel pointing out 
important properties of the solvers and their implementation 
in uG, as well as Anton Vadimovich Chizhov of the Physical- 
technical Institute loffe (Russian Academy of Sciences) in St. 
Petersburg, Russia for providing data for validation of the Borg 
Graham Channel model (Graham, 1999) used in conjunction 
with our hybrid method. 

SUPPLEMENTARY MATERIAL 

The Supplementary Material for this article can be found online 
at: http://www.frontiersin.org/journal/10.3389/fninf2014. 
00068/abstract 

Supplemental Movie 1 | Simulation of 200 ms real time, perspective 1. 
Supplemental Movie 2 | Simulation of 200 ms real time, perspective 2. 

REFERENCES 

AUbritton, N., Meyer, T., and Stryer, L. (1992). Range of messenger action of 
calcium ion and inositol 1,4,5-triphosphate. Science 258, 1812-1815. doi: 
10.1126/science.l465619 

Andrews, S. S., Addy, N. J., Brent, R., and Arkin, A. R (2010). Detailed simula- 
tions of cell biology with Smoldyn 2.1. PLoS Comput. Biol. 6:el000705. doi: 
10.1371/journal.pcbi.l000705 

Anwar, H., Hepburn, 1., Nedelescu, H., Chen, W., and De Schutter, E. (2013). 
Stochastic calcium mechanisms cause dendritic calcium spike variability. /. 
Neurosci. 33, 15848-15867. doi: 10.1523/JNEUROSC1.1722-13.2013 

AreUano, ]. 1., Benavides-Piccione, R., DeFelipe, ]., and Yuste, R. (2007). 
Ultrastructure of dendritic spines: correlation between synaptic and spine 
morphologies. Front. Neurosci. 1, 131-143. doi: 10.3389/neuro.Ol. 1.1.010.2007 

Ascoli, G. A. (2006). Mobilizing the base of neuroscience data: the case of neuronal 
morphologies. Nat. Rev. Neurosci. 7, 318-324. doi: 10.1038/nrnl885 

Bading, H. (1998). Transcription-dependent neuronal plasticity: the nuclear cal- 
cium hypothesis. Eur. J. Biochem. 267, 5280-5283. 

Balls, G. T., and Baden, S. B. (2004). "A large scale monte carlo sim- 
ulator for cellular microphysiology," in Proceedings of the 18th 
International Parallel and Distributed Processing Symposium (San 
Diego, CA: CaUfornia University). doi: 10.1046/j. 1432-1327.2000. 
01565.x 

Bastian, R, Birken, K., lohannsen, K., Lang, S., Neuss, N., Rentz-Reichert, H., et al. 

(1997). UG a flexible software toolbox for solving partial differential equations. 

Comp. Vis. Sci. 1, 27-40. doi: 10.1007/s007910050003 
Bentley, 1. L. (1975). Multidimensional binary search trees used for associative 

searching. Commun. ACM 18, 509-517. doi: 10.1145/361002.361007 
Blackwell, K. T., Wallace, L. J., Kim, B., Oliveira, R. E, and Koh, W. (2013). Modeling 

spatial aspects of intracellular dopamine signaling. Methods Mol. Biol 964, 

61-75. doi: 10.1007/978-l-62703-251-3_5 
Bower, J., and Beeman, D. (1997). The Book of GENESIS: Exploring Realistic Neural 

Models with the GEneral NEural Simulation System. New York, NY: Springer. 
Brandi, M., Brocke, E., Talukdar, H. A., Hanke, M., Bhalla, U. S., Kotaleski, 

J. H., et al. (2011). Connecting MOOSE and NeuroRD through MUSIC: 

towards a communication framework for multi-scale modeling. BMC Neurosci. 

12(Suppl. 1):P77. doi: 10.1186/1471-2202-12-sl-P77 
Broser, R 1., Schulte, R., Roth, A., Helmchen, E, Waters, Lang, S., et al. (2004). 

Nonlinear aniostropic diffusion filtering of three-dimensional image data from 

2-photon microscopy/. Biomed. Opt. 9, 1253-1264. doi: 10.1117/1.1806832 



Frontiers in Neuroinformatics 



www.frontlersin.org 



July 2014 I Volume 8 | Article 68 | 13 



Grein et al. 



1 D-3D hybrid modeling 



Burette, A. C, Lesperance, T., Crum, J., Martone, M., Volkmann, N., Ellisman, 
M. H., et al. (2012). Electron tomographic analysis of synaptic ultrastructure. 
/. Comp. Neurol 520, 2697-2711. doi: 10.1002/cne.23067 

Cannon, R. C, O Donnel, C, and Nolan, M. R (2010). Stochastic ion channel gat- 
ing in dendritic neurons: morphology dependence and probabilistic synaptic 
activation of dendritic spikes. PLoS Comp. Biol. 6:el000886. doi: 10.1371/jour- 
nal.pcbi.1000886 

Carnevale, N. T., and Mines, M. L. (2006). The NEURON Book. Cambridge UK, 

Cambrdige University Press, doi: 10.1017/CBO9780511541612 
Chen, X., Winters, C, Azzam, R., Li, X., Galbraith, ). A., Leapman, R. D., et al. 

(2008). Organization of the core structure of the postsynaptic density. Proc. 

Natl. Acad. Sci. U.S.A. 133,4453-4458. doi: 10.1073/pnas.0800897105 
Clapham, D. E. (2007). Calcium signaling. Cell 131, 1047-1058. doi: 

10.1016/j.cell.2007.11.028 
Cormen, T. H., Leiserson, C. E., and Rivest, R. L. (2011). Introduction to Algorithms. 

Cambridge, MA: MIT Press; Massachusetts Institute of Technology. 
Deuflhard, R, and Weiser, M. (2011). Numerische Mathematik 3. Berlin: De Gruyter 

Lehrbuch2011. 

Djurfeldt, M., Hjorth, J., Eppler, J. M., Dudani, N., Helias, M., Potjans, T. C, et al. 

(2010). Run-time interoperability between neural network simulators based on 

the MUSIC fi-amework. Neurinform 8, 43-60. doi: 10.1007/sl2021-010-9064-z 
Egelman, D. M., and Montague, R R. (1999). Calcium dynamics in the EC space 

of mammalian neural tissue. Biophys. J. 76, 1856-1867. doi: 10.1016/S0006- 

3495(99)77345-5 

Eggermann, E., Bucurenciu, I., Goswami, S. P., and Jonas, P. (2011). Nanodomain 
coupling between Calcium channels and sensors of exocytosis at fast mam- 
malian synapses. Nat. Rev. Neurosci. 13, 7-21. doi: 10.1038/nrn3125 

Eymard, R., Gallouet, T. R., and Herin, R. (2000). "Parabolic equations. Chapter 4," 
in The Finite Volume Method in Handbook Of Numerical Analysis, Vol. 7. eds P. 
G. Ciarlet and J. L. Lions (Elsevier), 713-1020. 

Gewaltig, M. O., and Diesmann, M. (2007). Nest (neural simirlation tool). 
Scholarpedial, 1430. 

Graham, L. (1999). Interpretations of Data and Mechanisms for Hippocampal 
Pyramidal Cell Models. New York, NY: Plenum Publishing Corporation. 

GriUo, A., Logashenko, D., Stichel, S., and Wittum, G. (2010). Simulation 
of density-driven flow in fractured porous media. Adv. Water Resour. 33, 
1494-1507. doi: 10.1016/j.advwatres.2010.08.004 

Hackbusch, W. (1985). Multi-Grid Methods. Leipzig: Springer. 

Hackbusch, W. (2003). Multi-Grid Methods and Applications in Springer 
Computational Mathematics. Leizpig: Springer. 

Hansen, S., Henning, A., NAd'gel, A., Heisig, M., Wittum, G., Neimiann, D., 
et al. (2008). In-silico model of skin penetration based on experimentally 
determined input parameters. Part i: Experimental determination of parti- 
tion and diffusion coefficients. Eur. J. Pharm. Biopharm. 68, 352-367. doi: 
10.1016/j.ejpb.2007.05.012 

Hayashi, Y, and Majewska, A. (2005). Dendritic spine geometry: functional 
implication and regulation. Neuron 46, 529-532. doi: 10.10I6/j.neuron.2005. 
05.006 

Helmchen, P., Imoto, K., and Sakmann, B. (1996). Ca2+ buffering and action 
potential-evoked ca2+ signaling in dendrites of pyramidal neurons. Biophys. J. 
70, 1069-1081. doi: 10.1016/30006-3495(96)79653-4 

Heppner, I., Lampe, M., Nagel, A., Reiter, S., Rupp, M., Vogel, A., et al. (2013). 
"Software framework uG4: parallelmultigrid on the hermit supercomputer," in 
High Performance Computing in Science and Engineering, eds W. E. Nagel, D. H. 
Kroner, and M. M. Resch (Berlin-Heidelberg: Springer), 434-449. 

Hines, M. L., and Carnevale, N. T. (2003). "The NEURON simulation environ- 
ment," in The Handbook of Brain Theory and Neural Networks, ed M. A. Arbib 
Vol. 2, (Cambridge MA: MIT Press), 769-773. 

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane 
currents and its appUcation to conduction and excitation in nerve. /. Physiol. 
117, 550-544. 

lerusalimschy, R., Celes, W., and de Figueiredo, L. H. (2013). Lua: The Programming 

Language. Rio: MIT License, 
lungblut, D., Queisser, G., and Wittum, G. (2011). Inertia based filtering of 

high resolution images using a gpu cluster. Comp. Vis. Sci. 14, 181-186. doi: 

10.1007/S00791-012-0171-2 
Kanamori, T. (2013). Compartmentalized calcium transients trigger dendrite prun- 
ing in DrosophUa sensory neurons. Science 340, 1475-1478. doi: 10.1126/sci- 
ence.1234879 



Katona, G., Kaszas, A., Turi, G. F., Hajos, N., Tamas, G., Vizi, E. S., et al. 
(2011). Roller coaster scanning reveals spontaneous triggering of dendritic 
spikes in CAl interneurons. Proc. Natl. Acad. Sci. U.S.A. 108, 2148-2153. doi: 
10.1073/pnas.l009270108 

Kerr, R. A., Bartol, T. M., Kaminsky, B., Dittrich, M., Chang, J. C, Baden, S. B., et al. 
(2008). Fast monte carlo simulation methods for biological reaction-diffusion 
system in solution and on surfaces. SIAM J. Sci. Comput. 30, 3127-3249. doi: 
10.1137/070692017 

Kits, K. S., Dreijer, A. M. C, Lodder, 1. C, Borgdorff A., and Wadman, W. J. 
(1997). High intracellular calcium levels during and after electrical discharges in 
molluscan peptidergic neurons. Neuroscience 79, 275-284. doi: 10.1016/S0306- 
4522(96)00651-3 

Korkotian, E., and Segal, M. (1999). Release of calcium from stores alters the mor- 
phology of dendritic spines in cultured hippocampal neurons. Proc. Natl. Acad. 
Sci. U.S.A. 96, 12213-12215. doi: 10.1073/pnas.96.21. 12068 

Lee, D. T. (1997). Worst-case analysis for region and partial region searches in 
multidimensional binary search trees and balanced quad trees. Acta Inf. 9, 
23-29. 

Maravall, M., Mainen, Z. R, Sabatini, B. L., and Svoboda, K. (2000). 
Estimating intracellular calcium concentrations and buffering without wave- 
length ratioing. Biophys. J. 78, 2655-2667. doi: 10.1016/S0006-3495(00) 
76809-3 

McDougal, R. A., Hines, M. L., and Lytton, W. W. (2013a). Reaction- 
diffusion in the NEURON simulator. Front. Neuroinform. 7:28. doi: 
10.3389/fnin£2013.00028 

McDougal, R. A., Hines, M. L., and Lytton, W. W. (2013b). Water-tight mem- 
branes from neuronal morphology files. /. Neurosci. Methods 220, 167-178. doi: 
10.1016/j.jneumeth.2013.09.011 

Migliore, M., Morse, T. M., Davison, A. P., Marenco, L., Shepherd, G. M., 
and Hines, M. L. (2003). ModelDB: making models publicly accessible 
to support computational neuroscience. Neuroinformatics 1, 135-139. doi: 
10.1385/NI:1:1:135 

Milner, B., L.R., S., and E.R., K. (1998). Cognitive neuroscience and the study of 

memory Neuron 20, 445-468. doi: 10.1016/30896-6273(00)80987-3 
Mount, D., and Arya, S. (2012). ANN Library 

Muha, I., Grillo, A., Heisig, M. Schoneberg, M., Linke, B., and Wittum, G. (2012). 
Mathematical modeling of process liquid flow and acetoclastic methanogenesis 
under mesophilic conditions in a two-phase biogas reactor. Biores. Technol. 106, 
1-9. doi: 10.1016/j.biortech.2011.11.087 

Muha, I., Nagel, A., Stichel, S., Grillo, A., Heisig, M., and Wittum, G. (2011). 
Effective diffusivity in membranes with tetrakaidekahedral cells and implica- 
tions for the permeabiUty of human stratimi corneum. /. Memb. Sci. 368, 18-25. 
doi: 10.1016/i.memsci.2010.10.020 

MuUer, D., Nikonenko, I., Jourdain, R, and Alberi, S. (2002). LTP, 
memory and structural plasticity. Cwrr. Mot Med. 2, 605-611. doi: 
1 0.2 1 74/ 1 566524023362041 

Murase, S., and Schuman, E. M. (1999). The role of cell adhesion molecules 
in synaptic plasticity and memory. Curr. Opin. Cell Biol. 11, 549—553. doi: 
10.1016/S0955-0674(99)00019-8 

Nagel, A., Hansen, S., Neumann, D., Lehr, C. M., Schaefer, U. F., Wittum, G., et al. 
(2008). In-silico model of skin penetration based on experimentally determined 
input parameters. Part ii: Mathematical modelling of in-vitro diffusion experi- 
ments, identification of critical input parameters. Eur. J. Pharma. Biopharm. 68, 
368-379. doi: 10.1016/j.ejpb.2008. 11.009 

Nagel, A., Heisig, M., and Wittum, G. (2009). A comparison of two- and three- 
dimensional models lor the simulation of the permeability of human stra- 
tum corneum. Eur. J. Pharm. Biopharm. 72, 332-338. doi: 10.1016/j.ejpb. 
2008.11.009 

Nagerl, U. V., Mody, I., and Vergara, J. L. (2000). Binding kinetics of calbindrn- 
d(28k) determined by flash photolysis of caged ca(2-H). Biophys. J. 79, 
3009-3018. doi: 10.1016/50006-3495(00)76537-4 

Oliveira, R., Terrin, A., Di Benedetto, G., Cannon, R. C, Koh, W., Kim, 
M., et al. (2010). The role of type 4 phosphodiesterases in generating 
microdomains of cAMP: large scale stochastic simulations. PLoS ONE 5:el 1725. 
doi: 10.1371/iournal.pone.0011725 

Popov, V. I., Kleschevnikov, A. M., Klimenko, O. A., Stewart, M. G., and Belichenko, 
p. V. (2011). Three-dimensional synaptic ultrastructure in the dentate gyrus and 
hippocampal area ca3 in the ts65dn mouse model of down syndrome. /. Comp. 
Neurol. 519, 1338-1354. doi: 10.1002/cne.22573 



Frontiers in Neuroinformatics 



www.frontiersin.org 



July 2014 I Volume 8 | Article 68 | 14 



Grein et al. 



1 D-3D hybrid modeling 



Pumplin, D. W., Reese, T. S., and Llinas, R. (1981). Are the presynaptic membrane 
particles the calcium channels? Proc. Natl. Acad. Sri. U.S.A. 11, 7210-7213. doi: 
10.1073/pnas.78.11.7210 

Queisser, G., Wittmarm, M., Bading, H., and Wittum, G. (2008). Filtering, recon- 
struction, and measurement of the geometry of nuclei from hippocampal 
neurons based on confocal microscopy data. /. Biomed. Opt. 13:014009. doi: 
10.1117/1.2829773 

Ray, S., and Bhalla, U. S. (2008). PyMOOSE: interoperable scripting in Python for 

MOOSE. Front. Neuroinf. 2:6. doi: 10.3389/neuro.ll.006.2008 
Resasco, D. C, Gao, R, Morgan, R, Novak, I. L., Schaff, J. C., and Slepchenko, B. M. 

(2012). Virtual cell: computational tools for modeling in cell biology. Wiley 

Interdiscip. Rev. Syst Biol. Med. 4, 129-140. doi: 10.1002/wsbm.l65 
Shewchuk, J. R. (2002). "What is a good linear finite element? - interpola- 
tion, conditioning, anisotropy and quality measures," Proceedings of the 11th 

International Meshing Roundtahle (Berkeley, OA). 
Si, H. (2009). TetGen: a quality tetrahedral mesh generator and three-dimensional 

delaunay triangulator. Int. J. Num. Methods Eng. 75, 856-880. 
Silver, R. A., Lamb, A. G., and Bolsover, S. R. (1990). Calcium hotspots caused by 

1-channel clustering promote morphological changes in neuronal growth cones. 

Nature 343, 751-754. doi: 10.1038/343751a0 
Spacek, J., and Harris, K. M. (1997). Three-dimensional organization of smooth 

endoplasmic reticulum in hippocampal CAl dendrites and dendritic spines of 

the immature and mature rat. /. Neurosci. 24, 4233^241. 
Tada, X, and Sheng, M. (2006). Molecular mechanisms of dendritic spine morpho- 
genesis. Curr. Opin. Neurohiol. 16, 95-101. doi: 10.1016/j.conb.2005.12.001 
Tai, C, Kim, S., and Schuman, E. (2008). Cadherins and synaptic plasticity. Curr. 

Opin. Cell Biol 5, 567-575. doi: 10.1016/j.ceb.2008.06.003 
Takahashi, A., Camacho, P., Lechleiter, D. J., and Herman, B. (1999). Measurement 

of intracellular calcium. Physiol. Rev. 79, 1089-1125. 
Thompson, ). R, Soni, B. K., and Weatherill, N. R (2012). Handbook of Crid 

Generation. Leipzig: CRC Press. 
Van Aelst, L., and Cline, H. (2004). Rho GTPases and activity-dependent 

dendrite development. Curr. Opin. Neurobiol. 14, 297—304. doi: 

10.1016/j.conb.2004.05.012 
Vogel, A., Reiter, S., Rupp, M., Nagel, A., and Wittum, G. (in press). uG 4 - A 

novel flexible software system for the simulation of PDE-based models on high 

performance computers. Comput. Vis. Sci. 
Wald, I., and Hvran, V. (2006). "On building fast kd-trees for ray tracing and 

on doing that in o(NlogN)," IEEE Symposium on Interactive Ray Tracing 

(Amsterdam), doi: 10.1109/RT.2006.280216 



Weinan, E., Engquist, B., and Huang, Z. (2003). Heterogenous multiscale method: 

a general methology for multiscale modeling. Phys. Rev. B 67, 367—450. doi: 

10.1 103/PhysRevB.67.092101 
West, A., Griffith, E., and Greenberg, M. (2002). Regulation of transcription 

factors by neuronal activity. Nat. Rev. Neurosci. 3, 921-931. doi: 10.1038/ 

nrn987 

Wikipedia (2012). K-d tree — Wikipedia, The Free Encyclopedia. 

Wils, S., and De Shutter, E. (2009). STEPS: modeling and simulationg com- 
plex reaction-diffusion systems with python. Front Neuroinf 3:15. doi: 
10.3389/neuro.ll.015.2009 

Wittmann, M., Queisser, G., Eder, A., Wiegert, J., Bengtson, C., Hellwig, A., 
et al. (2009). Synaptic activity induces dramatic changes in the geometry of 
the cell nucleus: interplay between nuclear structure, histone H3 phospho- 
rylation, and nuclear calcium signaling. /. Neurosci. 29, 14687—14700. doi: 
10.1523/INEUROSCI.l 160-09.2009 

Wolf, S., Grein, S., and Queisser, G. (2013). Employing NeuGen 2.0 to auto- 
matically generate realistic morphologies of hippocampal neurons and neu- 
ral networks in 3D. Neuroinformatics 11, 137-148. doi: 10.1007/sl2021-012- 
9170-1 

Xylouris, K. (2013). Generalized cable equation for signal processing in neurons. 
PhD Thesis. 

Xylouris, K., Queisser, G., and Wittum, G. (2007). A three-dimensional mathe- 
matical model of active signal processing in axons. Comput. Visual. Sci. 13, 
409-^18. doi: 10.1007/s00791-01 1-0155-7 

Conflict of Interest Statement: The authors declare that the research was con- 
ducted in the absence of any commercial or financial relationships that could be 
construed as a potential conflict of interest. 

Received: 09 January 2014; accepted: 28 June 2014; published online: 29 July 2014. 
Citation: Grein S, Stepniewski M, Reiter S, Knodel MM and Queisser G (2014) 1D-3D 
hybrid modeling — from multi-compartment models to full resolution models in space 
and time. Front Neuroinform. 8:68. doi: 10.3389/fninf.2014.00068 
This article was submitted to the journal Frontiers in Neuroinformatics. 
Copyright © 2014 Grein, Stepniewski, Reiter, Knodel and Queisser. This is an open- 
access article distributed under the terms of the Creative Commons Attribution License 
(CC BY). The use, distribution or reproduction in other forums is permitted, provided 
the original author(s) or licensor are credited and that the original publication in this 
journal is cited, in accordance with accepted academic practice. No use, distribution or 
reproduction is permitted which does not comply with these terms. 



Frontiers in Neuroinfoimatics 



www.frontiersin.org 



July 2014 I Volume 8 | Article 68 | 15 



